#include<stdio.h>
#include<stdlib.h>
#include<sys/time.h>
#include<math.h>
#include<string.h>

// Erros
#define ERROR_INVALID_FILE -1
#define ERROR_INVALID_MATRIX -2
#define ERROR_INVALID_VECTOR_SIZE -3


// Retorna o maior valor em modulo
float abs_max(float *vector, int size) {
    // Indice do maior elemento
    int abs_max_index = 0;
    // Iterador
    int i;

    // Tamanho invalido
    if(size <= 0) return ERROR_INVALID_VECTOR_SIZE;

    // Percorre o vetor
    for(i = 0; i < size; i++) {
        if(fabs(vector[i]) > fabs(vector[abs_max_index])) abs_max_index = i;
    }

    // Retorno
    return fabs(vector[abs_max_index]);
}

// Rotina principal

int main(void) {
    // Tempo de execucao
    struct timeval tvInitial, tvFinal;
    double tInitial, tFinal;

    // DEBUG
    int debug = 0;

    // Variaveis
    int J_ORDER, J_ROW_TEST, J_ITE_MAX;
    float J_ERROR;
    float **MA, *MB;
    float **MADivDiag, *MBDivDiag;
    float *Xn, *XnNext;
    float *Xdiff;
    float resultNext;
    char fileName[500];
    char filePath[500];
    FILE *fileInput;
	FILE *fileResult;
	FILE *fileDebug;

    // Iteradores
    int i, j;

    // Decremento
    float *decrement;

    // Erro inicial
    float errorNext;

    // Numero de iteracoes
    int numIterations = 0;

    // PASSO 1: ABERTURA DO ARQUIVO
    scanf("%s", fileName);

    // Tempo de execucao inicial
    gettimeofday(&tvInitial, NULL);
    tInitial = (double)(tvInitial.tv_sec) + (double)(tvInitial.tv_usec)/ 1000000.00;

	// Abertura do arquivo de entrada
    strcpy(filePath, "\0");
    strcat(filePath, "data/");
    strcat(filePath, fileName);
    fileInput = fopen(filePath, "r");
    if(fileInput == NULL) {
        printf("ERRO! Arquivo invalido!\n");
        return ERROR_INVALID_FILE;
    }

	// Abertura do arquivo de saida
    strcpy(filePath, "\0");
    strcat(filePath, "result_");
    strcat(filePath, fileName);
    fileResult = fopen(filePath, "w");
    if(fileResult == NULL) {
        printf("ERRO! Arquivo de saida invalido (Permissao negada)!\n");
        return ERROR_INVALID_FILE;
    }

	// Abertura do arquivo de debug
	if(debug) {
		strcpy(filePath, "\0");
		strcat(filePath, "debug_");
		strcat(filePath, fileName);
		fileDebug = fopen(filePath, "w");
		if(fileDebug == NULL) {
		    printf("ERRO! Arquivo de saida invalido (Permissao negada)!\n");
		    return ERROR_INVALID_FILE;
		}
	}

    //
    // PASSO 2: LEITURA DOS DADOS
    //

    // Ordem da matriz
    fscanf(fileInput, "%d", &J_ORDER);

    // Fila para ser avaliada
    fscanf(fileInput, "%d", &J_ROW_TEST);

    // Erro permitido
    fscanf(fileInput, "%f", &J_ERROR);

    // Numero maximo de iteracoes
    fscanf(fileInput, "%d", &J_ITE_MAX);

    // ALOCACAO DAS MATRIZES

    // Xn
    Xn = (float*) calloc(J_ORDER, sizeof (float));
    // Xn+1
    XnNext = (float*) calloc(J_ORDER, sizeof (float));
    // Xdiff
    Xdiff = (float*) calloc(J_ORDER, sizeof (float));

    // MA
    MA = (float**) malloc(J_ORDER * sizeof (float*));
    for (i = 0; i < J_ORDER; i++) {
        MA[i] = (float*) calloc(J_ORDER, sizeof (float));
    }

    // MB
    MB = (float*) calloc(J_ORDER, sizeof (float));

    // MADivDiag
    MADivDiag = (float**) calloc(J_ORDER, sizeof (float*));
    for (i = 0; i < J_ORDER; i++) {
        MADivDiag[i] = (float*) calloc(J_ORDER, sizeof (float));
    }

    // MBDivDiag
    MBDivDiag = (float*) calloc(J_ORDER, sizeof (float));

    // Decremento
    decrement = (float*) calloc(J_ORDER, sizeof (float));

    // LEITURA DAS MATRIZES

    // Matriz A (i x j)
    for (i = 0; i < J_ORDER; i++) {
        for (j = 0; j < J_ORDER; j++) {
            fscanf(fileInput, "%f ", &MA[i][j]);
        }
    }

    // Matriz B (1 x j)
    for (j = 0; j < J_ORDER; j++) {
        fscanf(fileInput, "%f", &MB[j]);
    }

    // DEBUG: Imprime os dados
    if(debug) {
        fprintf(fileDebug, "DADOS DO PROBLEMA\n");
        fprintf(fileDebug, "Ordem da matriz: %d\n", J_ORDER);
        fprintf(fileDebug, "Fila para ser avaliada: %d\n", J_ROW_TEST);
        fprintf(fileDebug, "Erro permitido: %f\n", J_ERROR);
        fprintf(fileDebug, "Numero maximo de iteracoes: %d\n\n", J_ITE_MAX);
        fprintf(fileDebug, "MATRIZ A\n");
        for (i = 0; i < J_ORDER; i++) {
            for (j = 0; j < J_ORDER; j++) {
                fprintf(fileDebug, "%f ", MA[i][j]);
            }
            fprintf(fileDebug, "\n");
        }
        fprintf(fileDebug, "\n");
        fprintf(fileDebug, "MATRIZ B\n");
        for (j = 0; j < J_ORDER; j++) {
            fprintf(fileDebug, "%f ", MB[j]);
        }
        fprintf(fileDebug, "\n\n");
    }

    //
    // PASSO 3: EXECUCAO DO ALGORITMO
    //

    // Encontra a matriz MADivDiag
    for (i = 0; i < J_ORDER; i++) {
        for (j = 0; j < J_ORDER; j++) {
            // Verificacao da diagonal
            if (MA[i][i] == 0) {
                printf("ERRO! Elemento da diagonal principal zerado!\n");
                return ERROR_INVALID_MATRIX;
            };
            // Diagonal principal
            if (i == j) MADivDiag[i][j] == 0;
            // Demais elementos
            else {
                MADivDiag[i][j] = MA[i][j] / MA[i][i];
            }
        }
    }

    // Encontra a matriz MBDivDiag
    for (i = 0; i < J_ORDER; i++) {
        // Verificacao da diagonal
        if (MA[i][i] == 0) {
            printf("ERRO! Elemento da diagonal principal zerado!\n");
            return ERROR_INVALID_MATRIX;
        }
        // Elementos
        MBDivDiag[i] = MB[i] / MA[i][i];
    }

    // DEBUG: Imprime os dados
    if(debug) {
        fprintf(fileDebug, "MATRIZ A*\n");
        for (i = 0; i < J_ORDER; i++) {
            for (j = 0; j < J_ORDER; j++) {
                fprintf(fileDebug, "%f ", MADivDiag[i][j]);
            }
            fprintf(fileDebug, "\n");
        }
        fprintf(fileDebug, "\n");
        fprintf(fileDebug, "MATRIZ B*\n");
        for (j = 0; j < J_ORDER; j++) {
            fprintf(fileDebug, "%f ", MBDivDiag[j]);
        }
        fprintf(fileDebug, "\n\n");
    }

    // Solucao inicial
    for(i = 0; i < J_ORDER; i++) {
        Xn[i] = MBDivDiag[i];
    }

    // Metodo iterativo
    do {
        // Calcula Xk+1[j]
        for(i = 0; i < J_ORDER; i++) {
            // Reseta o decremento
            decrement[i] = 0.0f;

            // Xk+1[i] = B*[i]
            XnNext[i] = MBDivDiag[i];

            // Calcula o decremento (sum(A*[i][j] * xk[j])
            for(j = 0; j < J_ORDER; j++) {
                decrement[i] += MADivDiag[i][j] * Xn[j];
            }

            // Xk+1[i] -= sum(A*[i][j] * xk[j]
            XnNext[i] -= decrement[i];

            // Xdiff[i] = | Xk+1[i] - Xk[i] |
            Xdiff[i] = fabs(XnNext[i] - Xn[i]);
        }

        // Calcula o valor aproximado resultado do processo
        resultNext = 0.0f;
        for(j = 0; j < J_ORDER; j++) {
            resultNext += MA[J_ROW_TEST][j] * XnNext[j];
        }

        // Calcula o erro
        errorNext = abs_max(Xdiff, J_ORDER) / abs_max(XnNext, J_ORDER);

        // DEBUG: Imprime os dados completos da iteracao
        if(debug) {
            fprintf(fileDebug, "ITERACAO %d\n", numIterations);
            fprintf(fileDebug, "RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resultNext, MB[J_ROW_TEST]);
            fprintf(fileDebug, "abs_max(Xdiff): %f\n", abs_max(Xdiff, J_ORDER));
            fprintf(fileDebug, "abs_max(XnNext): %f\n", abs_max(XnNext, J_ORDER));
            fprintf(fileDebug, "Error %.3lf\n\n", errorNext);
        }

        // Atualiza a solucao
        for(i = 0; i < J_ORDER; i++) {
            Xn[i] = XnNext[i];
        }

        // Atualiza o numero de iteracoes
        numIterations++;
    }
    while(errorNext >= J_ERROR && numIterations < J_ITE_MAX);

    // Imprime os resultados do processo no arquivo de resultado
    fprintf(fileResult, "Iterations: %d\n", numIterations - 1);
    fprintf(fileResult, "RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resultNext, MB[J_ROW_TEST]);

    // Tempo de execucao final
    gettimeofday(&tvFinal, NULL);
    tFinal = (double)(tvFinal.tv_sec) + (double)(tvFinal.tv_usec)/ 1000000.00;

    // Imprime o tempo de execucao no arquivo de resultado
    fprintf(fileResult, "Tempo de execucao: %.2lfs\n", (tFinal - tInitial));

	// Imprime os resultados do processo na tela
    printf("Iterations: %d\n", numIterations - 1);
    printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resultNext, MB[J_ROW_TEST]);

	// Encerra os arquivos de entrada e saida
	fclose(fileInput);
	fclose(fileResult);
	if(debug) fclose(fileDebug);

	// Desaloca os ponteiros alocados
	for(i = 0; i < J_ORDER; i++) free(MA[i]);
	free(MA);
	free(MB);
	for(i = 0; i < J_ORDER; i++) free(MADivDiag[i]);
	free(MADivDiag);
	free(MBDivDiag);
	free(Xn);
	free(XnNext);
	free(Xdiff);
    free(decrement);

    // Retorno
    return 0;
}
